#!/usr/bin/python
# GM_InterpRoutines.py

import math
import numpy


def linear(x1,y1,x2,y2):
    # Linear interpolation
    # Solved from y = m*x + b
    m = (y1 - y2)/(x1 - x2)
    b = y1 - m*x1
    return (m, b)


def quadratic(x1,y1,x2,y2,x3,y3):
    # Parabolic interpolation
    # Solved from y = a*x**2 + b*x + c
    a = ((x2-x1)*(y3-y1) + (y2-y1)*(x1-x3))/((x3**2-x1**2)*(x2-x1) + (x3-x1)*(x1**2-x2**2))
    b = (y2 - y1 + a*(x1**2 - x2**2))/(x2-x1)
    c = y1 - a*x1**2 - b*x1
    return (a, b, c)


def cubic(x1,y1,x2,y2,x3,y3,x4,y4):
    # Cubic interpolation
    # Solved from y = d*x**3 + e*x**2 + f*x + g
    # Way uglier than quadratic interpolation, due in part to the fact that I refused to solve by hand
    # d; simplify messy Mathematica solution for sanity
    x12 = x1**2*x2 - x1*x2**2
    x13 = x1**2*x3 - x1*x3**2
    x14 = x1**2*x4 - x1*x4**2
    x23 = x2**2*x3 - x2*x3**2
    x24 = x2**2*x4 - x2*x4**2
    x34 = x3**2*x4 - x3*x4**2
    dnum = -(y1*(x23 - x24 + x34) + y2*(-x13 + x14 - x34) + y3*(x12 - x14 + x24) + y4*(-x12 + x13 - x23))
    ddenom = (x1 - x3)*(x2 - x3)*(x2 - x4)*(x4 - x3)*(x1 - x2)*(x1 - x4)
    d = dnum/ddenom
    # e; exploit simplification & quadratic solution
    y1p = y1 - d*x1**3
    y2p = y2 - d*x2**3
    y3p = y3 - d*x3**3
    e = ((x2-x1)*(y3p-y1p) + (y2p-y1p)*(x1-x3))/((x3**2-x1**2)*(x2-x1) + (x3-x1)*(x1**2-x2**2))
    # f
    y1s = y1 - d*x1**3 - e*x1**2
    y2s = y2 - d*x2**3 - e*x2**2
    f = (y1s - y2s)/(x1 - x2)
    # g
    g = y1s - f*x1
    return (d,e,f,g)


def DeltasLimitsLinear(some1darray, clipends=True):
    '''Delimit the bins s.t. the limit is halfway between two "centers" provided by some1darray.

    If clipends == True, the limits of some1darray are the absolute limits and the two end bins are half as wide as they would otherwise be.
    
    A point of note.  Deltas must by definition be positive, but we don't know a priori whether some1darray runs from small to large numbers or the other way around, or whether it is all positive.  Likewise with Limits, some care must be exercised w.r.t. signs.  It IS assumed that some1darray is monotonic.

    Also, this is (secretly) approved for multi-dimensional inputs, not just 1d arrays, but it will always calculate the difference/gradient along the 0th axis.'''
    some1darray = numpy.array(some1darray, dtype=float)
    if type(clipends) != bool:
        print "Warning! Check your input type to GM_InterpRoutines.DeltasLimitsLinear. The optional input clipends should be boolean."
    # Bin widths
    Deltas = numpy.zeros(some1darray.shape)
    Deltas[1:-1] = 0.5*(some1darray[2:] - some1darray[:-2])
    if clipends == False:
        Deltas[0] = some1darray[1] - some1darray[0]
        Deltas[-1] = some1darray[-1] - some1darray[-2]
    if clipends == True:
        Deltas[0] = 0.5*(some1darray[1] - some1darray[0])
        Deltas[-1] = 0.5*(some1darray[-1] - some1darray[-2])
    Deltas = numpy.absolute(Deltas)
    # Bin limits
    shapetuple = some1darray.shape
    shapetuple = list(shapetuple)
    shapetuple.append(2)
    shapetuple = tuple(shapetuple)
    extradimensions = len(shapetuple) - 2
    dstring = ':,'*extradimensions
    # 
    Limits = numpy.zeros(shapetuple)
    #Limits[1:,%0] = 0.5*(some1darray[:-1] + some1darray[1:])
    #Limits[:-1,%1] = 0.5*(some1darray[:-1] + some1darray[1:])
    cmd = "Limits[1:,%s0] = 0.5*(some1darray[:-1] + some1darray[1:])" % dstring
    exec cmd
    cmd = "Limits[:-1,%s1] = 0.5*(some1darray[:-1] + some1darray[1:])" % dstring
    exec cmd
    if clipends == True:
        #Limits[0,0] = some1darray[0]
        #Limits[-1,1] = some1darray[-1]
        cmd = "Limits[0,%s0] = some1darray[0]" % dstring
        exec cmd
        cmd = "Limits[-1,%s1] = some1darray[-1]" % dstring
        exec cmd
    if clipends == False:
        #Limits[0,0] = some1darray[0] - 0.5*(some1darray[1] - some1darray[0])
        #Limits[-1,1] = some1darray[-1] + 0.5*(some1darray[-1] - some1darray[-2])
        cmd = "Limits[0,%s0] = some1darray[0] - 0.5*(some1darray[1] - some1darray[0])" % dstring
        exec cmd
        cmd = "Limits[-1,%s1] = some1darray[-1] + 0.5*(some1darray[-1] - some1darray[-2])" % dstring
        exec cmd
    return (Deltas, Limits)



def FillInGaps(some1darray, maxstep):
    '''A function to reduce the unevenness of sampling in some1darray.

    Leaves all existing points in place, just interjects new ones.

    Increasing and decreasing, positive and negative, are allowed; monotonic is assumed.'''
    maxstep = math.fabs(float(maxstep))
    some1dlist = list(some1darray)
    tiny = 1e-6
    newlist = [some1dlist[0]]
    for i in range(1, len(some1darray)):
        delta = some1dlist[i] - some1dlist[i-1]
        if math.fabs(delta) <= maxstep*(1 + tiny):
            newlist.append(some1dlist[i])
        else:
            # Round up, get number of evenly spaced partitions needed
            numberintervals = int(math.fabs(delta/maxstep) + 1.0 - tiny)
            intervalsize = delta/float(numberintervals)
            for j in range(0, numberintervals):
                newentry = some1dlist[i-1] + intervalsize*(j+1)
                newlist.append(newentry)
    newarray = numpy.array(newlist)
    return newarray


def NewEvenArray(old1darray, step, strictstep=False, offset=False):
    '''Creates a new array with the limits of old1darray and spacing step.
    
    If strictstep == True, step size is enforced even if it means overrunning the limits of old1darray, and those limits are overrun by the same amount at each end.  If strictstep == False, the limits of old1darray are enforced, and step may be shrunk in order to fit evenly into the available range.
    
    If offset == False, the limits are left where they are.  If offset == True, the values are offset so that bins centered on them with width step have upper and lower limits at the original limiting values.

    old1array is assumed monotonic.  However, whether it is increasing or decreasing, the output will be the same.'''
    tiny = 1e-6
    upperlimit = old1darray[-1]
    lowerlimit = old1darray[0]
    delta = upperlimit - lowerlimit
    if type(strictstep) != bool:
        print "The type of parameter 'strictstep' must be boolean."
    if type(offset) != bool:
        print "The type of parameter 'offset' must be boolean."
    if strictstep == True:
        # Round up
        numberofsteps = int(math.fabs(delta/step) + 1.0 - tiny)
        intervalsize = math.copysign(step, delta) # returns step w/sign of delta
        finaldelta = numberofsteps*intervalsize
        midpoint = 0.5*(upperlimit + lowerlimit)
        finalupper = midpoint + 0.5*finaldelta
        finallower = midpoint - 0.5*finaldelta
    if strictstep == False:
        # Round up
        numberofsteps = int(math.fabs(delta/step) + 1.0 - tiny)
        intervalsize = delta/float(numberofsteps)
        finalupper = upperlimit
        finallower = lowerlimit
    if offset == True:
        finallower = finallower + 0.5*intervalsize
        finalupper = finalupper - 0.5*intervalsize
        numberofsteps = numberofsteps - 1
    if offset == False:
        pass
    newarray = numpy.zeros(numberofsteps+1)
    for i in range(0,numberofsteps+1):
        newarray[i] = finallower + intervalsize*i
    return newarray


def OldToNewLinear(oldx, oldy, newx):
    '''Linearly interpolates between values of oldy at positions oldx to find values newy at positions newx.'''
    # Establish a "small" length for these arrays
    oldspaces = oldx[1:] - oldx[:-1]
    oldspaces = numpy.abs(oldspaces)
    oldspaces = oldspaces.min()
    newspaces = newx[1:] - newx[:-1]
    newspaces = numpy.abs(newspaces)
    newspaces = newspaces.min()
    tiny = 1.e-3*min(newspaces, oldspaces)
    # Fill in the new array
    shapetuple = list(oldy.shape)
    shapetuple[0] = len(newx)
    shapetuple = tuple(shapetuple)
    newy = numpy.zeros(shapetuple)
    for i in range(0, len(newx)):
        for j in range(0, len(oldx)):
            # If it is (nearly) equal to an old entry
            if math.fabs(newx[i] - oldx[j]) < tiny:
            #if numpy.less(numpy.abs(newx[i] - oldx[j]), tiny):
                newy[i] = oldy[j]
            # If it is before the beginning
            elif (j == 0) and (newx[i] < oldx[0]):
                x1 = oldx[0]
                y1 = oldy[0]
                x2 = oldx[1]
                y2 = oldy[1]
                (m, b) = linear(x1, y1, x2, y2)
                newy[i] = m*newx[i] + b
            elif j == (len(oldx)-1):
                # If it is beyond the far end
                if newx[i] > oldx[j]:
                    x1 = oldx[j-1]
                    y1 = oldy[j-1]
                    x2 = oldx[j]
                    y2 = oldy[j]
                    (m, b) = linear(x1, y1, x2, y2)
                    newy[i] = m*newx[i] + b
                # Prevents next condition from choking
                break
            # If it is really close to the next old entry, wait
            elif math.fabs(newx[i] - oldx[j+1]) < tiny:
                continue
            # If it is between two old entries
            elif (oldx[j] < newx[i]) and (oldx[j+1] > newx[i]):
                x1 = oldx[j]
                y1 = oldy[j]
                x2 = oldx[j+1]
                y2 = oldy[j+1]
                (m, b) = linear(x1, y1, x2, y2)
                newy[i] = m*newx[i] + b
    return newy



def DownSampleLinear(newlims, oldlims, oldstuff):
    '''Resamples oldstuff from oldlims to newlims along outer (0th) axis.

    Best that newlims and oldlims cover the same range - I haven't really accounted for alternatives.

    Assumes oldstuff entries have been multiplied by bin width.
    
    oldlims and newlims must have shapes (n1, 2) and (n2, 2), where n1, n2 may be any integer.  n1 must also be the length of oldstuff along outer axis.'''
    newlims = numpy.array(newlims)
    oldlims = numpy.array(oldlims)
    oldstuff = numpy.array(oldstuff)
    shapelist = list(oldstuff[0].shape)
    shapelist.insert(0, len(newlims))
    shapetuple = tuple(shapelist)
    newstuff = numpy.zeros(shapetuple, dtype=float)
    for i in range(0, len(newlims)):
        hilim = newlims[i,1]
        lolim = newlims[i,0]
        newwidth = hilim - lolim
        tempbuffer = numpy.zeros(oldstuff[0].shape, dtype=float)
        for j in range(0, len(oldlims)):
            oldhi = oldlims[j,1]
            oldlo = oldlims[j,0]
            oldwidth = oldhi - oldlo
            if oldhi < lolim:
                continue
            elif oldlo > hilim:
                continue
            # If old falls inside new
            elif (hilim >= oldhi) and (lolim <= oldlo):
                inside = 1.
                outside = 1.
                tempbuffer = tempbuffer + oldstuff[j]
            # If old encompasses new
            elif (hilim < oldhi) and (lolim > oldlo):
                outside = oldhi - oldlo
                inside = hilim - lolim
                tempbuffer = tempbuffer + oldstuff[j]*inside/outside
            # If old overlaps new top
            elif (hilim < oldhi) and (lolim <= oldlo):
                outside = oldhi - oldlo
                inside = hilim - oldlo
                tempbuffer = tempbuffer + oldstuff[j]*inside/outside
            # If old overlaps new bottom
            elif (hilim >= oldhi) and (lolim < oldlo):
                outside = oldhi - oldlo
                inside = oldhi - lolim
                tempbuffer = tempbuffer + oldstuff[j]*inside/outside
            
#            tempbuffer = tempbuffer + oldstuff[j]*inside/outside
        # And what if boundaries don't match up?
        # Well, screw that.  I rule that they will match up.
        newstuff[i] = numpy.array(tempbuffer)
    return newstuff
